﻿#include "changWeiFenSolver.h"

typedef __int64 MyInt;

/*
** y'= f(t,y)
** y(t0)=y0
** function f(double,double) is f in upper equ
** method is 4th level 4th order classical Runge-Kutta method
*/

bool changWeiFenSolver::SolveAndGetY_RungeKutta_4Order(double x0, double y0, double xMax, double maxStep_dh, std::function<double(double t, double y)> f, std::vector<CoordPoint>& vecResult, std::string& errorMsg)
{
	if ((xMax - x0) <= 0)
	{
		errorMsg = "error: xMax must be above x0.";
		return false;
	}
	MyInt numOfPoints = (xMax - x0) / maxStep_dh;
	if (numOfPoints == 0)
	{
		errorMsg = "error:xMax must >= (x0 + maxStep_dh).";
		return false;
	}

	if (maxStep_dh <= 0)
	{
		errorMsg = "error: maxStep_dh must >0";
		return false;
	}

	//x0 x1
	bool bJustRight = ((xMax - x0 - numOfPoints * maxStep_dh) / maxStep_dh) < 1e-13;
	if(!bJustRight)
	{
		numOfPoints += 1;
	}
	numOfPoints += 1;//加上初值点
	vecResult.resize(numOfPoints);

	vecResult[0].x = x0;
	vecResult[0].y = y0;

	for (MyInt i = 1; i < numOfPoints; i++)
	{
		if (i == numOfPoints - 1)//last point
		{
			vecResult[i].x = xMax;
		}
		else
		{
			vecResult[i].x = vecResult[i - 1].x + maxStep_dh;
		}
	}
	double dh = maxStep_dh;
	double k1 = 0;
	double k2 = 0;
	double k3 = 0;
	double k4 = 0;

	for (MyInt i = 1; i < numOfPoints; i++)
	{
		if (i == numOfPoints - 1)
		{
			dh = vecResult[i].x - vecResult[i - 1].x;
		}

		k1 = f(vecResult[i - 1].x, vecResult[i - 1].y);
		k2 = f(vecResult[i - 1].x + 0.5 * dh, vecResult[i - 1].y + 0.5 * dh * k1);
		k3 = f(vecResult[i - 1].x + 0.5 * dh, vecResult[i - 1].y + 0.5 * dh * k2);
		k4 = f(vecResult[i - 1].x + dh, vecResult[i - 1].y + dh * k3);

		vecResult[i].y = vecResult[i - 1].y + dh * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
	}
}